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Abstract 

The main contribution of this work is to construct higher than second order accurate 
total variation diminishing (TVD) schemes which can preserve high accuracy at non-sonic 
extrema with out induced local oscillations. It is done in the framework of local maximum 
principle (LMP) and non-conservative formulation. The representative uniformly second 
order accurate schemes are converted in to their non-conservative form using the ratio of 
consecutive gradient. These resulting schemes are analyzed for their non-linear LMP/TVD 
stability bounds using the local maximum principle. Based on the bounds, second order 
accurate hybrid numerical schemes are constructed using a shock detector. Numerical results 
are presented to show that such hybrid schemes yield TVD approximation with second or 
higher order convergence rate for smooth solution with extrema. 
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1 Introduction 

We consider the ID scalar conservation law associated to the conserved variable u{x,t), 

2tt(x, t) + t)) = 0, (x, p e M X M+ 

u{x, 0) = uo{x) 


( 1 ) 


where f{u) is a non-linear flux function. The numerical approximation for the solution of Q 
is done by the discretization of the spatial and temporal space into N equispaced cells Ij = 
\Xi_i,Xi,-i\, z = 0,1,... V of length Ax and M equispaced intervals n = 0,l,...,M 

of length At respectively. Let Xj = zAx and t"" = nAt denote the cell center of cell Ij and the 
time level respectively then a conservative numerical approximation for ([^ can be defined 

by 

A=v (2) 

where tt” = ri(xj,t"') and is the numerical flux function defined at the cell interface x-^i 


at time level n. The characteristics speed a{u) = associated with (j^ can be approximated 
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as 


Tpn rpn 

^i+1 



^i+1 “i 


du 


if / nf, 

if <+!=< 


( 3 ) 


where Fi = f{uf). In general due to non-linearity of Q, beyond a small finite time, even 
for a smooth initial data the evolution of discontinuities in the solution u{x, t) is inevitable. 
Therefore, it is required to have a conservative approximation of the solution with high accuracy 
and crisp resolution of such discontinuities with out numerical oscillations. Contrary to this 
need, most classical high order schemes despite of being linearly Von-Neumann stable give 
oscillatory approximation for discontinuities even for the trivial case of transport equation i.e., 
f{u) = au, 0 ^ a € R. Such oscillatory approximation can not be considered as admissible 
solution since it violets the following global maximum principle satished by the physically correct 
solution u{x,t) of Q i.e.. 


min(uo(x)) < u{x,t) < inax{uo{x)){x, t) G M x M"*". 


(4) 


In order to overcome these undesired numerical instabilities, various notion of non-linear stability 
are developed in the light of maximum principle Q . Examples of Maximum principle satisfying 
schemes are monotone schemes |42l H]. total variation diminishing (TVD) schemes |13[ [38l I45l 
[Ml 0 Ea ESI ST]. Some uniformly high order maximum-principle satisfying and positivity 
preserving schemes are Ensansj. There are other non-oscillatory schemes which do not strictly 
follow maximum principle but practically give excellent numerical results e.g.. Essentially non- 
oscillatory (ENO) and weighted ENO schemes see |35) and references therein. It is known 
that among global maximum principle satisfying schemes, the monotone and total variation 
diminishing (TVD) schemes experience difficulties at data extrema. On the one hand, such 
high order schemes locally degenerate to hrst order accuracy at non-sonic data extrema and on 
the other hand, even such a uniformly hrst order accurate schemes may exhibit induced local 
oscillations at data extrema. In this work the focus is on the construction of improved TVD 
schemes at smooth data extrema. 


2 Global maximum principle and data extrema 

The above global maximum principle Q satisfying monotone and TVD schemes have been of 
great interest mainly due to excellent convergence proofs for entropy solution |33| (22] and |281 03] 
respectively. The key idea is, any maximum principle satisfying scheme produce a bounded 
solution sequence and convergence follows due to compactness of solution sequence space |23|. 
It can be shown that monotone stable scheme TVD scheme ^ monotonicity preserving scheme 
(or Local extremum diminishing (LED)) scheme |11[ 1^. Unfortunately, monotone as well TVD 
schemes experience difficulty at data extrema. The monotone stability relies on monotone data 
and therefore a monotone scheme preserves the monotonicity of a data set by mapping it to 
a new monotone data set but fails to preserve the non-monotone solution region i.e., at data 
extrema. These monotone schemes are criticized mainly due to barrier theorem which state that 
a ’linear’ three point monotone scheme can be at most first order accurate [9]. Later, second 
order ’non-linear’ conservative monotone schemes are constructed using limiters but again by 
compromising on second order accuracy at extrema, e.g. |42| . The TVD stability mimics the 
maximum principle as it relies on the condition that global extremum values of solution must 
remain be bounded by global extremum values of initial solution. In !ni, Harten gave the 
concept of total variation diminishing scheme by measuring the variation of the grid values as 
follows 


2 





Definition 2.1. Conservative seheme 0 is called total variation diminishing if 


E 


-<+‘l < 


E 


-u,- 


( 5 ) 


where A+ttj = = ttj+i — Ui. 


Note that that the definition 2.1 is global as it is defined on the whole computational domain 
and ensures that global maxima or minima of initial solution uo{x) will not increase or decrease 
respectively. Such conservative TVD schemes are heavily criticized because, even if they are 
higher order accurate in most solution region, they give up second order of accuracy at non- 
sonic critical values of the solution |291 1280 We emphasize that these depressing results on 
degeneracy of accuracy of TVD method are given for conservative schemes and in the above 
global sense. More precisely the global nature of TVD definition (|^ allows shift in indices 
technique in ^ sign and is extensively used in different terms of the infinite sums in the TVD 
proofs of various schemes and results in the literature including the following one due to Harten 
HI]. 

Lemma 2.2. A conservative scheme in Incremental form (I-form) 

■”+‘ = < + a.+ . W+i - <) - /3.:_. W - <_,) 


u- 


( 6 ) 


is TVD iff a -,1 > 0, /3j,i > 0 and ot-.i + /3j, i < 1, Vi. fSBf 


2.1 Degenerate accuracy at extrema: 


In |28| . proof for degeneracy to first order accuracy at non-sonic critical points of solution i.e., 
points u*{x,t) s.t. / {u*(x,t)) ^ 0 = u’^{x,t) is mainly based on modified equation analysis 
and a conservative semi-discrete version of Lemma 2.2 In [29] . using a trade off between second 


order accuracy and TVD requirement along with shift in indices technique, it is shown that 
second order accuracy must be given up by a conservative TVD scheme at non sonic critical 
values Ui = u*{xi,t) which corresponds to extreme values i.e., [u{xi + Ax,t) — u{xi,t)].[u{xi,t) — 
u{xi — Ax, t)] < 0 ^ f (ui). It is also worthy to note that problem of degenerate accuracy by 
modern high resolution TVD schemes is also due to their construction procedure. For example, 
the numerical flux function of flux limiters or slope limiters based TVD schemes is essentially 
design in such a way that it reduces to first order accuracy at extrema and high gradient region 
by forcing limiter function to be zero see IMllE] and references therein. This makes it impossible 
for a limiter based TVD schemes to achieve higher than first order accuracy at solution extrema 
as well at steep gradient region |23]. Thus every high order TVD (in global sense ^ ) scheme 
suffers from clipping error and cause flatten approximation for smooth extrema though they 
sharply capture discontinuities [20] . 


2.2 Induced local oscillations: 

Apart from compromise in uniform high accuracy, it is notable that global maximum principle 
satisfying monotone and TVD schemes do not necessarily ensure preservation of non-monotone 
data set i.e., for a data set with extrema as demonstrated in Figure [^b). In particular first 
order monotone and TVD schemes with large coefficient of numerical viscosity can allow 
the occurrence of induced oscillations at data extrema and formation of new local extremum 
values as shown in Figure |^a). This phenomena of generation of local oscillations at extrema is 
reported and analyzed for well known monotone and TVD three point Lax-Friedrichs scheme in 
[H |26l |25] similar to Figure [^a). It is interesting to note that two point monotone and TVD 

^Sanders also defined the total variation by measuring the variation of the reconstructed polynomials and such 
TVD schemes can be uniformly high order accurate IHIIZI- 
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(a) (b) 

Figure 1: Induced oscillations may occur for non-monotone data extrema, (a) Monotone 
stability rely on monotone data sequence, (b) For non-monotone data set TVD definition Q is 
satisfied though updated approximation is oscillatory. 




(a) 


(b) 


Figure 2: Numerical approximation of Linear transport equation (35) with impulsive initial 
data, (a) By monotone and TVD Lax-Friedrichs scheme: Induced oscillations (b) First order 
two point upwind scheme: Absence of induced oscillations. 


upwind scheme does not exhibit induced oscillations in Figure j^b). 

The main aim of this work is to construct uniformly non-oscillatory shock capturing mono¬ 
tone and TVD methods with high accuracy for non-sonic smooth extremsj^ In order to achieve 
it, a non-conservative formulation is done using framework of local maximum principle (LMP) 
with the help of gradient ratio parameter. The rest of paper is organized as: For completeness, 
in section]^ local maximum principle (LMP) stability is defined for two points schemes. It is 
shown in Lemma 3.2 that away from sonic point LMP stability implies global monotone and 
TVD stability. In section]^ we analyze representative uniformly second order schemes in non- 
sonic region for their TVD (or LMP) stability by converting them into two point schemes. This 
yields computable bounds for the stability of these scheme and are presented as main results in 
Theorems 4.1 4.7 These obtained TVD bounds ensure for second rrder TV stable approximation 
for smooth solution with non-sonic critical point of extreme nature. In section]^ hybrid schemes 
are designed using the TVD bounds and a shock detector technique. Numerical results are given 
to support the theoretical results and claim. Conclusion and future work is discussed in section 

El 


^It shows improved TVD approximation in the region of degenerate accuracy reported in 
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3 Local maximum principle (LMP) stability 

It is clear from the above discussion that the global maximum principle Q satisfying monotone 
or TVD stability experience difficulty in the presence of non-monotone data with extrema. These 
two local phenomena i.e, induced oscillations and degeneracy in accuracy at non-sonic extrema 
by monotone and TVD schemes motivate us to look for their non-linear stability locally at non- 
sonic extrema. Consider the following local maximum principle (LMP|^for scalar conservation 
law 0, 

:,r) if f'{u) > 0, (7a) 

:,r) if f\u) < 0. (7b) 

In case of two point schemes, initial solution data will always be monotone as either u{xa,t) < 
u{xb,t), Xa 7^ Xb or vice verse thus the LMP condition Q reduces to 

min(nf_i,Mn < < max«_i,nf) if/(n) > 0, (8a) 

< max(n",n(Yi) if/(^) < 0- (8b) 


min u{x, 

e) 

< 

u{x,e+^) 

< 

max 

Xi—i<X<Xi 





Xi — i<X<Xi 

min u{x, 

e) 

< 

n(x,r+^) 

< 

max 

Xi<X<Xi+l 





Xi<X<X^J^l 


Thus away from sonic point u* i.e., a-_|_i = f (u*) ^ 0 define. 

Definition 3.1. A numerical scheme is LMP stable 0 if it can he written as 

if 0 < Aa” 1 < 1, 


= 


Cu"^ + 'Du^_^ 
Cu^ + Vu2+i 


if — 1 < Aa” 1 < 0, 


(9) 


where coefficients Cand P are real functions such that C>t),P>{),C + P = l. 

Scheme ([^ is essentially a convex combination of two point values of n(:, t”) thus ensures that 
the updated solution value of u{xi,t'^~^^) will remain be bounded by both point values without 
introducing of new local maxima-minima. Note that two point first order upwind scheme is a 
natural example of LMP stable scheme with coefficients 


C = 1 - Aa..!,!? = Ao._i if 0 < A/«) < 1 
C = 1 — = Ao^.i if — 1 < A/ (n”) < 0. 


(10a) 

(10b) 


This justifies the non-occurrence of local oscillation by first order upwind scheme in Figure |^b). 
From definition 3.1 it follows. 

Lemma 3.2. Local maximum principle stable scheme 0 is global total variation diminishing 
stable. 

Proof. Using relation C = 1 — P, rewrite ^ in the non-conservativ^ half incremental form 

( 11 ) 


as 


= 


<-P«-<_i) if 0<Aa-^, <1, 
< + if -l<Aa-/, <0, 


where from Definition 3.1 0 < C < 1 and 0 < 2? < 1. On appropriately choosing one of the 
coefficients a or /3 zero in I-form Q, the Lemma 2.2 shows approximation by half I-form 0 
is global TVD. □ 


^Also known as upwind range condition |2U| _ _ 

^Note that for problems with constant f'{u) e.g. linear transport equation (351, the form (111 is conservative. 


Also one can obtained a conservative approximation from 


§ 


by defining V suitably such as V = ^ ^ ^ q < 

A-u" 


\f'{u") < 1 as in |17) . In the persent work the coefficient T> comes out in such a form which results in to a 
non-conservative approximation. 
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Remark 1. The LMP stability is defined only in non-sonic region therefore by Lemma [3. 2 [ LMP 
stability implies TVD stability in non-sonic region i.e., away from sonic point / (u*) ^ 0. In 
this setting, from next section onward until stated the term LMP/TVD stability implies global 
TVD stability (|^ away from sonic point. 


4 Bounds on high order TVD accuracy 


In this section using definition 3.1, it is shown that second order total variation diminishing 


approximation is possible for the solution with non-sonic extreme critical points. It follows 
from the LMP stability bounds given for the representative second order accurate Lax-Wendroff 
(LxW), Beam-Warming (BW) and Fromm schemes respectively for scalar problem Q. Let 
the stencil of r -|- s -|- 1 point scheme locally does not contain sonic point u*{x,t) 

i.e., / (tt*) 7 ^ 0 and characteristics speed at local cell interfaces is non-zero. Note that the 
case of degenerate characteristic speed i.e., a-,i = 0 or/and a-i = 0 are not interesting as 
these schemes do not necessarily preserve their uniform order of accuracy. We also consider 
the wave speed slpit ai*" i -|- a~ i = a -, i such that ai*" , > 0 and a~ , <0. Note that for 

l-\-- 7.-1-.. 


> 0 


®+2 

1 > 0 , or, 1 =0 whereas a-^i < 0 




i+i 


^i+ 


1 = a^_^i < 0 , 0^1 = 0,. After 


- 2 --r 2 '2 '2 --r 2 2 

dropping the superscript for time level n, following function definitions and notations are used 

in the rest of the presentation. Define the smoothness parameter as 


f = r±(F,) = 


1 T \af. ) 






i 


lTAa-,)A±T: 


i± 


( 12 ) 


where the flux split F F^ = Fi is consistent with wave split and given by 

^+1 - (13) 


(t{x) = sgn{x) = 


(14) 


Here the superscript ± sign of rj denotes the positive/negative sign of wave speed. Also define 
the signum function 

+1 if X > 0 , 

— 1 if X < 0 . 

In order to analyze the local non-linear stability of considered schemes we choose practically 
viable CFL like condition 

0 < A max I/'(n) I < 1 (15) 


Note that the choice A max |/ {u)\ = 0 / (u) corresponds to the case of degenerate character- 

U 

istic speed or steady state case. 

4.1 Centered Lax-Wendroff scheme 


Theorem 4.1. Away from sonic point and under CFL condition (15), the second order accurate 


Lax- Wendroff scheme for scalar conservation law 0 is TVD in the solution region where 


rf G -cx), Ki Aa.^i 


*^2 


U ( 71 ( ) , CX) 


where numbers < 0, 71 > 0 depends on CFL number for linear stability given by 

1 — xa{x) 


KlIXl = - 


1 -V xa(x)) ’ 


(16) 
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(17) 


71(x) = 


xcr{x) 


2 + xa{x) 

Proof. Consider the numerical flux function of Lax-Wendroff (LxW) scheme 


1 1 


A+Ui. 


(18) 


In order to ensure non sonic region, let the characteristics speed is locally non-zero. Since LxW 
uses three point centred stencil [xi_s,Xj+r],r = s = 1 , it suffice to assume x a-_i > 0 . 


Case f'{u) > 0: Let a-^i > 0, then the conservative approximation using (18) can be 
written as 

Aa^^i / s Aa,_ 


= Ui- 


(^1 - Aa-^i^ A+Ui H-^ (^1 Aa-_i^ A^m 


which can be written in the following non-conservative half Incremental form 

A+Ui 


= Uj — 


Aa., 1 

'1-Aa 


i+i 


2 / A_rti 


+ 


Aa._i 


1 -|- A a- 1 
* 2 


A_Ui. 


From Lemma 3.2 half I-from (20) will be TVD if, 

A+Ui 


0 < 

which reduces to. 


Aa- 




- Aa._i (1 -I- Aa._i ) < Aa._,_i (1 - Aa 


A+rti 


*+5/ A-Ui - 

Note that Xa.iil — Aa,_i) > 0 under discrete CFL condition, 
* 2 * 2 

0 < A maxa^_|_i < 1, 


Hence inequality (22) can be written as 


or 


(1 + Aa._p ^ Afli+i) A+ai ^ 2 - Aa._i (1 + Aa._i) 

(1 - Aa._i) “ a-_i{l - X a-_i) A^m ~ A a._i (1 - A a._i) 


(1 + Aa._i) ^ ai+i(l-Aa.+ i) A+i(i ^ 2 + Ao._i 


Ao_i 
* 2 


(1 —Aa-i) a-i (1 — Aa-i) A_r(. 

‘2 2 2 

Using definition ^ and flux wave split ( |13| ), Inequality (25) becomes 

(1-l-Aa^i) (1 “ Aa+ 1 ) ^+_p+ 2-|-Aa^i 


“2 ^ _ ^+ 2 ' ^ 


(1-Aa+i) (1 - Aa+ 1 ) A^F+i Aa+ 

* 2 * 2 

Inequality (26) on inversion yields, 

r+ < 4(Aa+0 OR r+ > 7i+(Aa+ 0, 

2 2 

where 

(1 - Aa:^i)AlFj (1-Aa+i) 




A 


r- = 


(1-Aa+JA+Fi 


K~[ (Aa 


* 2 


(1 + Aa._i) 

* 2 


and 7 ^ (Aal^ 


(19) 


( 20 ) 


< 2 , 

( 21 ) 

(l + Aa.-i) 

( 22 ) 


(23) 


(24) 

(25) 

(26) 

( 27 ) 


2 +Aa+ 
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Case /(u) < 0: Let < 0, then the non-conservative I-form can be written as 


Ua — Ui 


Aa, 


A o- 1 


, \ A_nj 

^ V + a;^ 


A+Wj. 


Using Lemma 3.2 half I-from (28) will be TVD if. 

A a •, 1 / V A a. 


0 < -- 


i+i 


1 — A o, 


i+i 




A-Uj 

A+Wi 


< 1 


which can be written as 

AOi+l (l - A<..+ .) < -An,. . (l + Aa._ .) ^ < 2 + Aa,^. (l - Aa,^. 

The discrete CFL condition for a{u) < 0 is, 

— 1 < Xa-_^_i < 0, Vi. 


(28) 

(29) 

(30) 

(31) 


Therefore the qnantity — Aa-_(_i (1 -|-Aa-_(_i) > 0. Divide Ineqnality (30) by it and nsing (j^ 
yields, 


(l-Aa,^..) ^ (l + Aa._.) 


-W A_K 


1 + Xa-^1 


1 + AOi_,i I 


< - 


2 + Aa._^i (l - Aa.^i^ 


Aaj+i (1 + Aa._^i 


(32) 


or nsing flnx wave split (13) 




1 An. 1 


< 


1 A n. I 


A-Pr 


A a 1 — 2 


1 P Xu 


*+3 


A I U- 


< 


j+i 


A a. 


(33) 


i+i 


Ineqnality (33) on inversion yields. 


r. < Ki Aa 1 OR r. >71 Aa . . 


(34) 


where r ■ = 


1 ) ^+^i 


1 + Aa._i ) A_F- 


■> ( H+l ) = 


— ( 1 -|- Aa., 


1 — Aa 1 


and 71 ( Aa.^i ) = 


—Aa 


*+l 


2 — A a 


*+3 


Condition (27) and (34) completes the proof. 


□ 


4.1.1 LxW on Linear problem: Every extrema is non-sonic. 

In order to see the improvement in the TVD approximation at non-sonic extrema, consider the 
linear transport eqnation 


d d 

—a(a:, t) + a—{u{x, t)) = 0, a / 0 


(35) 


In this case the smoothness parameter (12) rednces to rf = 


— — -. Note that at point of 

A±Ui 


extrema the measnre of smoothness is negative i.e, for transport eqnation (35), every < 0 


implies a non-sonic extreme critical point. Following resnlt follows from Theorem 4.1 
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Corollary 4.2. Under the linear stability eondition 0 < A|a| < 1, the seeond order aceurate Lax- 
Wendroff seheme for (35) is total variation diminishing where ri G (oo, Ki(aA)) U (71 (aA), 00). 




Figure 3: (a) : LMP/TVD stability region (shaded) for Lax-Wendroff scheme (b): Zoomed view. 


In Figure]^ the behavior of CFL number aA on dependent parameters ki and 71 is shown. 
Note that when Xa —)• 0"*“ the parameter 71 —)• O"*" whilst when Xa —>■ 1“ the parameter ki —>-0“. 
In particular, for a|A| = 1, definitions (16), (17) yield non-TVD interval [ki = 0,71 = ^]. Note 
that under linear CFL condition 0 < Aa < 1, Ki(Aa) G (—1,0] and 7 i(Aa) G (0,1/3). The 
following result give CFL independent TVD bounds for LxW scheme [?]• 


Corollary 4.3. The Lax-Wendroff seheme for ( f35| ) is total variation diminishing under the 
linear stability eondition 0 < Ajaj < 1, if Vi =G (00, —1) U (^, 00) 

Thus it, can also be concluded from corollary |4.2| and shaded TVD region for LxW scheme 
in Figure 1^ that exeept for Vi G [^1,71] the seeond order aecurate LxW seheme yields TVD 
approximation for all solution region ineluding extreme points with < 0. More precisely. 


4.2 Upwind Beam-Warming scheme 


Theorem 4.4. Away from sonic point and under CFL eondition (15), seeond order accurate 
Beam-Warming scheme is TVD for scalar conservation law 0 in the solution region where 




( Aa± J ,72 ( Aa.^, 


where parameter 72 and K 2 defined as, 

K2{x) 


— (2 — xa{x)) 
xa{x) 


(36) 


72 (x) 


3 — xa{x) 
1 — xa(x) 


(37) 


Proof Case f'{u) > 0: In this case, BW stencil use grid points [xj_s,s = 2,r = 0 

thus to ensure locally non-sonic region it suffice to let a _k > 0,k = —1,1, 3. The numerical 

* 2 

flux of Beam-Warming scheme is. 


jpn,BW 

Vi 


i / \ 

T) H-^ (^1 - Aa-_1 j A.Uj, a-_^i > 0. 


(38) 
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Resulting non-conservative half I-form can be written as 


= Uj — 


Att-_r / \ Att-u 

" '3-Aa._i 1-^ (1-Aa„. 


A-Ui-\ 
5 / A-Uj 


A_ttj (39) 


For half I-from (39) to be TVD, from Lemma 3.2 


0 < 


An„-_i/ \ Att- 3 

^ '3-Aa._i)-^ (1-Aa._ 


A_nj_i 

3 I —^-< 1. 


* 2 / A_m 


(40) 


Compound Inequality (40) can be written as 


Aa-i (3 — Aa-i) — 2 < Aa-_3(1 — Aa-s)——^ < Aa,_i (3 — Aa,_i) (41) 
*2' *2" *2' *2' A_n,' 2" * 2" 


Under CFL condition (23), i.e., 0 < Aa^_|_i < 1, Vi quantity Aa^_i(l — Xa-_i) is positive, 


hence (40) can be written as. 


Aa._i(3 - Ao._i) - 2 ^ (1 - Aa-.s)^ (3-Aa._i) 

Aa-_i (1 - Aa-_i) “ Oj.i (1 - Aa._i) A^m “ (1 - Aa._i) 


or 


Aa._i - 2 ^ A_Ui-i 

“ a-_i (1 - Aa-_i) A_nj “(l-Aa-_i) 


Aa 1 
^ 2 


On using flux wave split (13) we get 


1 - 2 ^ (1 - Aa._3) A_F)ti ^ 

Aa +1 “ (r-AaA 1 ) A_F)+ “ (1 - Aa+,) 


Which, using (12) can be written as. 


K 2 ( Aa+ 1 1 < r+ 1 < 72 ( Aa+ i 


(42) 


where K 2 and 72 are dehned in (36) and (37). 


Case / (u) < 0: In case of negative characteristics speed, BW use stencil [xj_s, Xj+r]) s = 
0, r = 2] thus to ensure locally non-sonic region it suffice to let a-^k > 0,k = —1,1, 3. In 
case of negative wave speed the Beam-Warming flux is given by. 


j^n^BW ^ + Aa.+i) A+Ui+i, a.+ i > 0. 

Resulting non-conservative half I-form is. 


=Ui + 


A a, 


*+f / 1 , X 


A+Ui+i 


3 + Aa._^i 


(43) 


A+Uj. (44) 


Condition for (44) to be TVD is 


0< Aa.+i (l + Afli+l) ^^-Aa,+ i (3 + Aa,+ i) <2 


A+Ui 


(45) 
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A a,-j_ 1 ( 3 + A a,-j_ 1 ) < A s (1 + A a^_|_ s ) —^^— < 2 + A a^_^_ i (3 + A o^_|_ i ) (46) 


A+Ui 




^i+i 


Note under CFL condition (31) —1 < Aa^_|_i < 0, Vi, Aa^_|_i ^1 + Aa^_|_i^ is negative. 
Compound Inequality (46) reouced to 


2 + Aa,+ i (3 + Aa,^i) ^ a^+3 (l + Aa^^a) A+m+i ^ (g + Aa^+i) 






A+Ui 


1 + Aa._^i 


(47) 


or 


on using (13) 


+ ^ ^ «i+i (l + Aa.+i) ^ (3 + AQi+i) 


A o, 


i+| a._^i (l + Aa._^i^ (l + Aa._^i^ 


A a 1 +2 
^+2 

Xa~ 1 


< 


1 + A a., 3 


A+^i 


(l + Ao.+ i) 




< 


3 + A a 


*+l 


1 A u 


*+l 


(48) 


(49) 


Using (12), (36) and (37) Inequality (49) can be written as 


>^2 


(50) 


□ 


4.2.1 BW on Linear problem 


Corollary 4.5. Beam-Warming scheme for (35) is total variation diminishing under the linear 
stability condition 0 < A|o| < 1, if the smoothness parameter G [K 2 (Aa), 72 (Aa)] 




(a) (b) 

Figure 4: (a) : LMP/TVD stability region (shaded) for Beam-Warming scheme (b): Zoomed 

view. 

In Figure TVD region for Beam-Warming scheme is shown as shaded region. It can be 
deduced that 
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1. Aa —)• O"*" Beam-Warming scheme yield TVD approximation for ^ (—00,3]. 

2. Xa —7- 1“ Beam-Warming scheme yield TVD approximation for ^ [“I;00). 

Note that for A|a| G (0,1] parameters K 2 G (—00,—1] and 72 G [3,00). Following CFL number 
independent weaker TVD bounds can be concluded 


Corollary 4.6. The Beam-Warming scheme for (35) is total variation diminishing under the 
linear stability condition 0 < A|o| < 1, if = [—1,3]. 


4.3 Fromm’s scheme 

A less ocsillatory and second order accurate scheme is obtained by using a simple average of 
LxW and BW flux i.e., 


pn,FROMM ^ 1 f pn,LxW , pn,BW 

i+| “ 2 V 


(51) 


From Theorem 4.1 and Theorem 4.4 following result can be proved, 


Theorem 4.7. Away from sonic point and CFL condition (15), the Fromm’s scheme corre¬ 
sponding to flux (51) for scalar conservation law 0 is TVD in the solution region where 




—00, Ki ( Aa* 1 
2 


and 




«2 ( Aa± 1 


U ( 71 ( ) , cx) 


,72 ( 


where parameter ki, 7 i,K 2 CLnd 72 are defined in and (p^ respectively. 


5 Hybrid high order LMP/TVD stable schemes 


It follows from the Theorem 4.1 4.4 and Theorem 4.7 that it is possible to achieve second or 
higher order TVD approximation for most solution region including non-sonic exterma where 
Tj < 0. In order to demonstrate it numerically, we construct hybrid schemes using a mono- 
tone/TVD scheme as complementary conservative scheme (CCS). The following hybrid 
schemes are the natural choice which satisfies the LMP/TVD bounds obtained in previous sec¬ 
tion and thus ensures a LMP/TVD approximation. The second order accurate LMP/TVD 
schemes use second order LxW and BW schemes in the region of their LMP/TVD stability 


using bounds on smoothness parameter in Theorem 4.1 and 4.4 respectively, otherwise use a 
conservative conservative scheme (CCS). 


5.1 Centered scheme: LW-CCS 

1 : if rf < Ki OR rf > 71 then 
2 : Update LxW scheme 

3: else 

4: Update ^ CCS. 

5: end if 


5.2 Upwind Scheme: BW-CCS 


if [rf^i > K 2 and < 72 ) then 
Update F- BW scheme 

else 

Update ^ CCS. 

end if 
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5.3 Centred-Upwind scheme: FLWBW-CCS approximation 


This LMP/TVD stable scheme can be obtained nsing Fromm scheme in its region of LMP/TVD 
stability nsing bonnds in |4.7| along with schemes 5.1 and|5.2|as follows 


if (rf < Ki OR rf > 71 ) AND > K 2 AND < 72 ) then 
Update Fromm’s scheme 

else if rf < OR rf > 71 then 

Update ■<— LxW scheme 

else if {rf^^ > K 2 AND rf--^ < 72 ) then 
Update •<— BW scheme 

else 

Update ^ CCS. 

end if 


Note that the scheme 5.1 5.3 are non-conservative as they are based on TVD conditions on the 
smoothness parameter dedced from the non-conservative form of stndied scheme^ Therefore 
they captnre the steady shock accnrately bnt may prodnce moving shock at wrong location see 
resnlts in Fignre|^a). Note that incorrect shock location by scheme [5.1| is legging whereas by 
scheme |5.2| it is leading to the exact shock location in Fignre|^a). It is interesting to see the 
scheme 5.3 cancels the leading and legging errors and gives exact shock location in Fignrejl^a). 
This phenomena of yielding wrong moving shock location by non-conservative schemes along 
with the shock correction criteria is well explained in |14| . The idea for shock correction is to 
apply locally a shock capturing conservative scheme in the vicinity of discontinuity using a shock 
detector. It is therefore, to capture the moving shock correctly, the following hybrid approach 
can be used. 


5.4 Shock Correction: SC-LW-CCS, SC-BW-CCS, SC-FLWBW-CCS hybrid 
schemes 

if Shock region then 

Update ^ use CCS, 

else 

Update •<— with either of algorithm 

end if 


5.1te.3 


5.5 Extension to system of hyperbolic conservation laws 

Consider the hyperbolic systems of conservation law in one dimensions, 

d <9 ^ , 

-u+^F(u)=0. 


(52) 


where u is vector of conserved quantities , i = 1,2,.. .1 and F is the vector flux function. The 
above proposed schemes for scalar case are extended to non-linear systems ( |52[ ) in the natural 
manner using flux vector splitting and average flux Jacobian matrix A = F (u) of the flux 
function. In particular for the numerical results presented in next section, the Steger-Warming 
flux vector splitting is used for ID and 2D systems. The average Jacobian matrix is computed 
as follows. 




u 


i+l 


+ U' 


®except for equations having constant characteristic speed e.g. linear transport problem. 
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In order to compute the TVD bounds, the n—characteristic speeds associated with system 52) 


o' 1 

*+3 


^ + l I 
7 7 

"4+1-'' 




-ul^O 


i) else, 




(53) 


where is the spectrum of eigen values of • In above computation (53) the nonphysical 


discrete wave speed caused by numerical overflow in case of ~ rtj are corrected using 

following way which is similar to the wave speed correction technique proposed in |15) . It is 
done as, 


— 


i+i 


- <^71 


if lo-^ 


i+i 


> CFr 


a. 


j _ 


— CTn 


if lo-^ 


• ill ^ ^mim 
*■’“2 


(54) 

(55) 


where amax and amin refer to the local maxima and minima of the magnitudes of characteristic 


speeds associated with system (52). For example, the one dimensional Euler equations has the 


eigenvalues u, u± c where u and c denotes fluid velocity and the speed of sound respectively. 
In this case, we dehne 


o'max = max (max (|m|, |ti — c|, |n + c|)j, max (|ri|, |n — c|, |ri + c|)i+i), (56) 

Urnin = max (min (|m|, |rt — c|, |tt + c|)i, min (|u|, |tt — c|, |rt + c|)j+i). (57) 


5.6 Shock Sensor 

In order to locate the presence of discontinuities, a shock detector proposed in m with some 
modihcation is used. A brief detail on the shock switch is given below for the sake of completeness 
of the discussion on numerical implementation. 


Step 1: Check multigrid ratio check 


MR{i,h) = 


Tc{i,h) 


Tpii, h) + e 

where Tc{i,2h) and Tp{i,h) are the (4*^j5*^ and 6*^) order truncation error sum on a 
coarse (with N/2 grid points) and hne grid (with N points) respectively. The derivatives 
in this step are calculated by by sixth order compact scheme proposed in |23) . The small 
parameter 0 < e << 1 is used to avoid division by zero. 

Step 2: Calculate the local ratio check at the grid point Xi which has multigrid ratio MR{i, h) < 4, 


LR{i) = 


‘■R 


- U 




where = 3?ri — 4rij_i +ttj _2 and = 3?ri — 4rij_|_i +ttj +2 left and right slope respectively 
at grid point Xj. 


Step 3: Use a cutoff value 6 G (0,1] to create a shock switch (SS) on the result of step 2. i.e., 

SS{i) = 


0, if LR{i) < S i.e., data is locally smooth around grid point Xi, 
1, if LR{i) > 5 i.e., data is discontinuous around grid point Xi. 


Note that the above shock detector has parameters e and 6 which governs the sensitivity of 
shock switch. It is observed in numerical computations that for larger value of parameters e.g., 
e = 1 X 10“^, 6 = 0.8 above shock switch is less sensitive for mild shock or sharp turns and detects 
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only strong shocks whereas e = 1 x 10“®,(5 = 0.2 detect corners and mild shock along with the 
strong shock. Also in case of non-linear systems, it is observed that small oscillations may arise 
in the vicinity of shock depending on the choice of shock parameters e, 6. It is therefore, to make 
it robnst and less prone to parameters €,6 a slight modihcation is done as follows. 

Step 4: Treat neighbouring grid point Xi±i in discontinuity region if SS{i) = 1 i.e., use SS{i±l) = 
1 if SS{i) = 1. 


6 Numerical Results 

In this section numerical results are presented for various benchmark scalar and system test prob¬ 
lems in both one and two dimensions. Different smooth as well discontinuous initial conditions 
are taken to show the performance of schemes in section in terms of accuracy and disconti¬ 
nuity capturing respectively. Numerical results show that the proposed hybrid scheme, due to 
improved accuracy at extrema and steep gradient region, nicely approximates the smooth region 
of solution with crisp resolution for rarefaction, contact and shock discontinuities. Moreover it 
produces the total variation diminishing numerical approximation. 


6.1 Linear transport equation: every extrema is non-sonic 

Consider the linear transport equation 


ut + Ux = 0, u{x, 0) = uo{x) 


(58) 


with periodic boundary condition. The exact solution of (58) equation convects with out chang¬ 


ing the initial shape of uo{x) and is given by u{x, t) = uo{x — t). Note that in general number of 
extrema are hnite and clipping error due to degenerate accuracy at smooth extrema by existing 
high order monotone and TVD method is visible only after a long period of time in form of 
approximation of smooth extrema with corners or flatten prohle see Figure |^a). Also due to 
degenerate accuracy at extrema and steep gradient region their erratic convergence rate can be 
seen only after a long time see Table It is therefore, probably the transport equation (58) 
the only test which can be used to check the large time performance of any method. Since the 


problem (58) is linear thus discontinuities present in the solution does not represents shock or 


rarefaction therefore scheme 5.3 can be directly applied for this test problem with out shock 


switch. The numerical computation for problem (58) is done by using the hrst order upwind 


scheme as complementary conservative scheme (CCS) in hybrid scheme 5.3 and results are shown 
by legend method. 


6.1.1 Accuracy check: smooth initial condition 


Consider 58) along with the the following three different initial conditions which comprises of 
smooth extrema, monotone region with mild as well sharp turn. 

i Smooth extrema 


u(x,0) =sin(7ra:), x G [—1, 1]. (59) 

The initial prohle consists smooth extreme points at x = which preservingly convects 
to the right direction. 


ii Smooth extrema with monotone data region 

u{x,0) = sin^(7rx), x G [0, 1]. (60) 


This initial condition is taken from m and has a smooth extrema at x = 0.5 with 
monotone solution regions with mild turn towards the bottom. 
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iii Smooth extrema with steep gradient region 


u(x,0 ) 


-1 

giz^ X G [—1 : 1] 

0 else 


(61) 


This initial condition has a smooth extrema at x = 0. Compared to initial condition (60) 
it has high gradient monotone region with sharp tnrns towards the bottom where r —)• 0+ 
or r >> 1 respectively. This is a good test to see the degenerate convergence rate for any 
limiter based TVD scheme. 


For smooth initial data nnmerical solntion plots are given at large time level Tg The conver¬ 
gence rate of scheme |5.4| is given at varions time instance with varying CFL nnmber to show 
the robnst and higher than second order convergence rate of the method |5.3[ In Fignrej^a), 


nnmerical solntion obtained by method corresponding to IC (59) is compared with high order 


TVD Lax-Wendroff flnx limited method (LxWFLM) |3T] with compressive Snperbee limiter 
In Fignre|^and Fignrej^ approximate solntion is given for transport problem corresponding to 
initial conditions ( [60] ) and ( |61[ ) respectively. The total variation of the compnted solntion ob¬ 
tained by method is compared with nniformly second order LxW, BW and Fromm schemes 
respectively for all three IC’s in Fignrej^b), |^b) and[^c). 

From the nnmerical resnlts in it is evident that problem of flattening of smooth ronnd shaped 
solntion prohle is removed dne to improved approximation of extreme points. Moreover it can 
be observed and Fignre [^a) that this improvement more visible as Aa —?■ 1 which snpport the 
improved TVD region for extrema of LxW as discnssed in the Corollary |4.2[ Total variation 


plots show that total variation for designed scheme 5.3 is decreasing whilst nniformly second 
order LxW and BW schemes do not prodnce TVD solntion thongh it remain total variation 
bonnded (TVB). 

In Table discrete maximnm error convergence rate is given for scheme [5(3] as error 
is the best indicator for checking the performance of any scheme in terms of clipping error dne 
to drop in accnracy at smooth extrema. In Table and Table error convergence rates are 
given in terms of and L°° error for method with different choice of CFL and time for linear 
test corresponding to initial conditions (60) and (61) respectively. The nnmerical resnlts show 
that the designed scheme 5.3 shows higher than second order convergence rate independently 
of the choice of CFL nnmber or hnal time T. Also dne to improved approximation of extrema 
and steep gradient region, the nsed method yields smooth approximation with ont clipping error 


which snpport the Corollary 4.2 and 4.3 


T=2 


T= 

=30 


CFL=0.5 

CFL=0.95 

CFL= 

0.5 

CFL=0.95 

N 

L°° error 

Rate 

L°° error 

Rate 

L°° error 

Rate 

L°° error 

Rate 

20 

9.1927e-03 


1.4026e-03 


6.3628e-02 


1.1716e-02 


40 

1.9456e-03 

2.240 

2.6640e-04 

2.396 

1.0892e-02 

2.546 

2.2095e-03 

2.407 

80 

3.7656e-04 

2.369 

5.6365e-05 

2.241 

1.8543e-03 

2.554 

4.4726e-04 

2.305 

160 

7.0744e-05 

2.412 

1.1332e-05 

2.314 

3.3400e-04 

2.473 

8.6187e-05 

2.376 

320 

1.2599e-05 

2.489 

2.6879e-06 

2.076 

4.1812e-05 

2.998 

1.4613e-05 

2.560 

640 

3.2485e-06 

1.955 

3.1057e-07 

3.113 

3.5142e-06 

3.573 

4.5598e-07 

5.002 


Table 1: Consistent higher than second order 
corresponding to initial condition (59). 


convergence rate with the mesh refinement 


®We consider large time Ts as the time level when corners are visible in the approximation of smooth extrema 
by high order TVD methods such as in [TTl [S] |38l |45l . 
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(b) at T = 2 with CFL = 0.95. 


Figure 5: Solution of (58) with IC (59) using N = 80: (a) Flatten approximation for smooth 
extrema by LxWFLM whereas proposed scheme 5.4 preserve solution with smooth extrema with 
out introducing corners, (b) Comparison of total variation. 




* 

LxW 


BW 

o 

Fromm 

X 

method 






time 


(a) 


(b) 


Figure 6: Solution for linear equation (35) with IC (60) (a) at time T = 20, 
Smooth approximation of extrema with reduced clipping error: (b) TV plot up to T 
N = 80, CFL = 0.95. 


N = 100, 
= 2 using 


6.1.2 Discontinuous initial condition 


In this test is taken from |6] originally used by Harten in |12) . Initial solution is complex in 
nature which contains parts of smooth solution, mix discontinuities, discontinuities of derivative 
in the interval [—1,1]. In Figure [^a), numerical results is given by proposed method and in 
Figure j^b) the total variation plot of the computed solution by method is given. 


' 2x + 2-sm(37r(x-0.5))/6 if -1 < x < 0.5, 
(0.5 — x)sm(1.57r(x — 0.5)^) if —0.5 < x < 1/6, 
|szn(27r(x — 0.5))| if l/6<x<=5/6, 

_ 2x-2-sm(37r(x-0.5))/6 if 5/6 < x <= 1. 


(62) 


It can be seen in Figure]^ that the proposed scheme [5^ yields TVD solution with crisp capturing 
of the discontinuities with out clipping error at smooth extrema. The error convergence rate is 
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CFL=0.5 


CFL= 

=0.95 


N 

Error 

Rate 

L°° error 

Rate 

Error 

Rate 

error 

Rate 

40 

7.6962e-02 


6.1348e-03 


1.3179e-02 


1.4315e-03 


80 

1.7350e-02 

2.149 

1.1039e-03 

2.474 

3.4723e-03 

1.924 

2.7228e-04 

2.394 

160 

3.8066e-03 

2.188 

1.9596e-04 

2.494 

8.7436e-04 

1.990 

5.2159e-05 

2.384 

320 

7.3723e-04 

2.368 

2.4368e-05 

3.008 

2.0210e-04 

2.113 

8.6669e-06 

2.589 

640 

1.1693e-04 

2.656 

4.3613e-07 

5.804 

3.2508e-05 

2.636 

1.1141e-07 

6.282 

1280 

2.8601e-05 

2.031 

5.4189e-08 

3.009 

8.0887e-06 

2.007 

1.3877e-08 

3.005 

2560 

7.0685e-06 

2.017 

6.6891e-09 

3.018 

2.0182e-06 

2.003 

1.7315e-09 

3.003 


Table 2: Order of convergence with the mesh refinement at T 
condition (60). 


20 corresponding to initial 




(a) (b) 


Figure 7: Solution for linear equation (35) with IC (61) at time T = 6 and N = 


0.95. (a) Smooth approximation of extrema with no clipping error by Scheme 5.1 
ison of total variation. 


100, CFL = 
(b) Compar- 


I-Order Upwind 

Second order TVD scheme 

N 

Error 

Rate 

L°° error 

Rate 

Error 

Rate 

error 

Rate 

80 

2.5248e-01 

0.583 

1.6091e-02 

1.482 

2.2319e-02 

1.314 

2.4725e-03 

2.294 

160 

1.6210e-01 

0.639 

5.8392e-03 

1.462 

1.6468e-02 

0.439 

1.3464e-03 

0.877 

320 

9.7673e-02 

0.731 

2.2646e-03 

1.367 

5.2704e-03 

1.644 

2.0531e-04 

2.713 

640 

5.6346e-02 

0.794 

8.2206e-04 

1.462 

2.3125e-03 

1.188 

8.7812e-05 

1.225 

1280 

3.1355e-02 

0.846 

2.8328e-04 

1.537 

1.0757e-03 

1.104 

4.1735e-05 

1.073 

2560 

1.6889e-02 

0.893 

9.2241e-05 

1.619 

4.1654e-04 

1.369 

1.0950e-05 

1.930 


Table 3: Order of convergence using I order upwind and LxW flux limited TVD method with 
Superbee limiter at T = 6,CFL = 0.5 corresponding to initial condition (61). 


not shown for this discontinuous test problem as discontinuities can only be approximated with 
at most first order of accuracy |36) 
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CFL=0.5 


CFL= 

=0.95 


N 

Error 

Rate 

L°° error 

Rate 

Error 

Rate 

L°° error 

Rate 

80 

4.3378e-02 

1.654 

4.9870e-03 

2.139 

1.1681e-02 

1.646 

2.2055e-03 

2.239 

160 

1.2992e-02 

1.739 

1.3187e-03 

1.919 

3.2801e-03 

1.832 

5.0305e-04 

2.132 

320 

3.2156e-03 

2.014 

2.6333e-04 

2.324 

8.1011e-04 

2.018 

8.1325e-05 

2.629 

640 

6.1867e-04 

2.378 

3.8218e-05 

2.785 

1.7801e-04 

2.186 

1.3173e-05 

2.626 

1280 

1.1353e-04 

2.446 

4.1315e-06 

3.210 

4.5829e-05 

1.958 

1.8927e-06 

2.799 

2560 

1.7115e-05 

2.730 

3.6114e-07 

3.516 

1.2376e-05 

1.889 

2.5213e-07 

2.908 


Table 4: Consistent higher than second order of convergence with the mesh refinement atT = 6 
corresponding to initial condition (61) 




(a) 


(b) 


Figure 8: (a) High resolution oscillation free solution for test problem 62 by scheme 5.4 for 


data CFL = 0.95, N = 160, T = 2.0. (b) Total variation decreasing plot of method 


6.2 Non-linear case: Burgers equation 

Consider the Burgers equation 


ut + 



= 0, —a < X < b 


(63) 


with initial condition uq{x) and periodic boundary conditions. It is the non-linear nature of the 


equation (63) that even for smooth initial condition, the solution of (63) eventually develops 


discontinuities like rarefaction and shocks after breaking time Tj, given by 


n = 


-1 


min3;(uo(a:)) 


(64) 


Also the unique sonic point for Burgers equation (63) is n* = 0. It is therefore, Burgers equation 


is a good test to check the performance of any scheme for smooth as well discontinuous solutions 
profile at pre and post-shock time Th respectively. In the numerical computation FORCE scheme 


is used as CCS in schemes 5.1 to |5.3| and their respective hybrid shock corrected analog in scheme 

[Ql 
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6.2.1 Shock correction moving shock 


We first consider the following discontinuous initial condition to show the non-conservative 
nature of schemes 5.1[ 5.2|and 5.3 as they yield solution with wrong location of moving shock. 


u{x, 0) = 


1 , 

0 , 






else. 


(65) 


The solution corresponding to IC (65), develops a rarefaction wave and a moving shock which 
corresponds to initial discontinuities at x = —1/3 and x = 1/3 respectively. In Figure [^a) re¬ 
sults obtained by second order Lax-Wendroff (LW-FORCE) and Beam-Warming (BW-FORCE) 


schemes 5.1 and 5.2 respectively are given. It is clear that shock location given by LW-FORCE 


is legging behind whereas by BW-FORCE it is leading ahead of exact shock location. Results 
by shock corrected schemes SC-LW-FORCE and SC-BW-FORCE as described in 5.4 are given 
in l^b) which show correct location of shock is recovered with out loosing crisp resolution of 
left rarefaction. Numerical results in Figure [Tw a), obtained by scheme 5.3 (FLWBW-FORCE) 


shows that shock is crisply captured at right location with high resolution for bottom and top 
of left rarefaction. In Figure [T^b), results by SC-FLWBW-FORCE are given which show little 
dissipative resolution for shock which is due to approximation by dissipative FORCE scheme in 
the vicinity of shock. 



-1 -0.5 0 o.s 1 



(a) 



-1 -o.s 0 0.5 1 



(b) 


Figure 9: Solution at T = 0.8 using CFL = 0.8, N = 80, (a) Wrong location of moving shock 
using Scheme 5.1 and 5.2 (b) Shock correction using shock switch. 


6.2.2 Accuracy Check: Smooth Initial conditions 

Consider three different smooth initial conditions (IC) along with corresponding breaking time 

n 

rt(x, 0) = 0.1 -h sm'‘(7rx), x G [0, 1], Ti, = 0.27803225 (66) 

?x(x, 0) = 0.1 -k exp-"^^ X G [-2 : 3], = 0.65669683. (67) 
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(a) (b) 


Figure 10: Solution at T = 0.8 using CFL = 0.8, N = 80, (a) Scheme 5.3 give correct 
shock location high resolution for left rarefaction and crisp capturing for the moving shock, (b) 


Solution by Scheme 5.4 


1 4 

u{x, 0) = -(1 + sm(7rx)), x G [—1 : 1], = —. 

4 TT 


( 68 ) 


IC (66) and (67) does not contain any sonic point whereas IC ( |68[ ) has a sonic point at x = —0.5 
since n(—0.5,0) = 0. The solution corresponding to IC (66) and (67) develop a moving shock 
followed by a rarefaction fan whereas the moving shock corresponding to IC (68) is separated by 
two rarefaction fans. In Figurej^ and [T^ the pre and post-shock solution of Burgers equation 
obtained by the shock corrected hybrid Scheme [5.4| (SC-FLWBW-FORCE) corresponding to IC 
(66), (67) and (68) respectively are given. The total variation plots are also given for different 
choices of CFL number Amax„|/(n)|. In Table to Table and errors are shown 

at pre-shock time using different CFL number. Results show that the hybrid scheme nicely 
approximated pre-shock solution with out clipping error and does not introduce induced oscilla¬ 
tions near shock in the post-shock solution. Moreover purposed method yields a total variation 
diminishing solution and shows a consistent convergence rate between second and third order in 
both the norms. 




(a) 


(b) 


Figure 11: Burgers equation solution using SC-FLWBW-FORCE corresponding to IC (66): (a) 


No clipping error for smooth extrema as well near shock zone CFL = 0.95, iV = 50, (b) Effect 
of CFL on total variation of computed solution 
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N 


CFL 

=0.6 



CFL= 

=0.9 



error 

Rate 

L°° error 

Rate 

Error 

Rate 

L°° error 

Rate 

10 

3.1074e-02 


1.2720e-02 


2.7767e-02 


1.0627e-02 


20 

1.3874e-02 

1.1634 

4.6798e-03 

1.4426 

5.8287e-03 

2.2521 

2.5638e-03 

2.0514 

40 

2.0530e-03 

2.7566 

3.9157e-04 

3.5791e 

1.0398e-03 

2.4868 

2.2157e-04 

3.5324 

80 

4.0690e-04 

2.3350 

6.9600e-05 

2.4921e 

2.3176e-04 

2.1657 

3.1460e-05 

2.8162 

160 

8.6630e-05 

2.2317 

1.2510e-05 

2.4760e 

5.2120e-05 

2.1527 

3.9000e-06 

3.0120 

320 

1.4700e-05 

2.5591 

1.2000e-06 

3.3820e 

1.0380e-05 

2.3280 

4.1000e-07 

3.2498 

640 

2.4800e-06 

2.5674 

5.0000e-08 

4.5850e 

2.4200e-06 

2.1007 

5.0000e-08 

3.0356 


Table 5: 


Convergence rate corresponding to IC (66) at time T = Ti,l2 


CFL=0.45 


CFL. 

=0.95 


N 

L^error 

Rate 

L°°error 

Rate 

L^error 

Rate 

error 

Rate 

10 

3.1120e-01 


1.0335e-01 


2.5190e-01 


9.7632e-02 


20 

6.4051e-02 

2.2805 

2.7574e-02 

1.9061 

4.7036e-02 

2.4210 

2.2831e-02 

2.0964 

40 

8.7373e-03 

2.8740 

2.9299e-03 

3.2344 

1.2098e-02 

1.9590 

5.2950e-03 

2.1083 

80 

1.9428e-03 

2.1690 

4.6987e-04 

2.6405 

2.3084e-03 

2.3899 

5.5543e-04 

3.2530 

160 

4.1491e-04 

2.2273 

6.1350e-05 

2.9371 

4.2048e-04 

2.4568 

5.5800e-05 

3.3153 

320 

9.0750e-05 

2.1928 

6.7500e-06 

3.1841 

8.6830e-05 

2.2758 

5.3900e-06 

3.3719 

640 

2.0870e-05 

2.1205 

7.3000e-07 

3.2089 

1.9470e-05 

2.1569 

5.6000e-07 

3.2668 


Table 6: Third order 
Tb = 0.65669683 


convergence rate corresponding to IC (67) at pre-Shock time Tb/2, 




Figure 12: Solution corresponding to IC (67): (a) No clipping error for smooth extrema as well 
near shock zone CFL = 0.8, N” = 80, (b) Effect of CFL on total variation diminishing plot of 
computed solution. 


6.3 Buckley Leverett Equation 

Consider Buckley-Leverett equation which has convex-concave flux. This equation physically 
represents the flow of a mixture of oil and water through a porous medium. 


du df{u) 
dt dx 


(69) 
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Figure 13: (a) Comparison of computed solution of corresponding to IC (68) at Tb = l/dvr 

CFL = 0.6, N = 80 (b) Total variation diminishing plot up to f = 2Tb 


CFL=0.45 

CFL=0.9 

N 

L1 ei-j-or 

Rate 

L°° error 

Rate 

error 

Rate 

L°° error 

Rate 

10 

3.9495e-02 

0.0000 

1.0467e-02 

0.0000 

1.9593e-02 

0.0000 

4.7957e-03 

0.0000 

20 

7.2108e-03 

2.4534 

1.8537e-03 

2.4974 

4.6493e-03 

2.0753 

1.1326e-03 

2.0821 

40 

2.1212e-03 

1.7653 

3.7326e-04 

2.3122 

1.0845e-03 

2.1000 

2.4774e-04 

2.1927 

80 

4.9995e-04 

2.0850 

7.0710e-05 

2.4002 

2.7441e-04 

1.9826 

4.2200e-05 

2.5535 

160 

9.2620e-05 

2.4324 

1.7250e-05 

2.0353 

6.6870e-05 

2.0369 

8.3800e-06 

2.3322 

320 

2.2230e-05 

2.0588 

4.1900e-06 

2.0416 

1.6580e-05 

2.0119 

1.7400e-06 

2.2679 

640 

5.5900e-06 

1.9916 

1.0400e-06 

2.0104 

4.1900e-06 

1.9844 

3.6000e-07 

2.2730 


Table 7: Convergence rate for test case corresponding to IC (68) at pre-shock time t = 2/7r 
CFL = 0.8 


The flux function is given by, 


f{u) = 


V? + a(l — uY 

Here a is viscosity ratio and u represents the saturation of water and lies between 0 and 1 

6.3.1 One moving shock 


(70) 


Consider equation (69) with a = h and initial condition 


u{x, 0) = 


1, X < 0, 


0, X > 0. 

The solution involves one single moving shock followed by an rarefaction wave. 

6.3.2 Two moving shock 


Consider equation (69) with a = 4 and subject to initial condition 


u{x, 0) = 


1, —0.5 < X < 0, 


(71) 


(72) 


0 , elsewhere. 

The solution involves two moving shocks, each followed by an rarefaction wave. In numerical 


simulation flux limited high resolution LxW TVD scheme m is used in hybrid scheme |5.4| as 
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CCS. The results corresponding IC ( |71[ ) and (72) are given in Figure[l^a) and[l^b) respectively. 
Results show that the proposed scheme sharply captures both the fast and slow shocks. The 
rarefaction waves are also approximated with high resolution. 




Figure 14: Numerical solution using CFL = 0.8 (a) N = 80 at time t = 0.75 (b) N = 100 at 
time t = 0.4; Sharp resolution for rarefaction fans and slow and fast moving shocks. 


6.4 ID Euler Equation 

The ID Euler equations of the Gas dynamics is given by 

( P \ ( \ 

where u = I pu 1 and F(u) = I pu^ +p 1 denotes vector of conservative variables and 

\ E ) \ {E+p)u ) 

conservative fluxes respectively. Variables p, u and p represents density, velocity and pressure 
respectively . The total energy e is defined by, 


e = 


7-1 2 


(74) 


where 7 is the ratio of specihc heat coefficients. We consider the four shock tube problems mod¬ 
eled by (73) to check the robustness of proposed scheme in section 5.5 These shock tube tests 


check any method in capturing the contact and shock discontinuity along with non-oscillatory 
high resolution approximation for smooth extrema. In all the numerical test a simple high 
resolution TVD flux limited centered (FLIC) HOI E] scheme with MINBEE limiter is used as 
CCS in |5.4| We denote results by this scheme by FLWBW-FLIC instead SC-FLWBW-FLIC. 


Numerical results are compared with FLIC to see the improvement in capturing the solution 
prohle by FLWBW-FLIC. Note that, the MINBEE limiter satisfies the universal TVD stability 
region given in |8j and therefore robustly works for both positive and negative characteristics 
speed associated with system (73) 
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Density (kg/m'^) 


(75) 


6.4.1 Shu-Osher shock tube test [3j 

, (3.857143,2.629369,10.3333) x <-4.0, 

[P,u,p)- (i + o.2sin(5x),0,1) x >-4.0. 

This test depicts shock interaction with a sine wave in density. The main challenge in this case 
is to capture both the complex small-scale smooth flow and shocks. In Figure results are 
presented an compared with FLIC scheme. It is evident from zoomed hgure mb) that the 
FLWBW-FLIC yields oscillation free approximation for shock with higher resolution compared 
to FLIC for complex oscillatory solution region about [0.5, 2.5]. It also capture the smooth region 
in around [—3, 0.5] with out clipping or flattening error which is due to improved approximation 
of smooth extrema and steep gradient region. 


Plot of Density vs Position Piot of Density vs Position 




(a) 


(b) 


Figure 15: Numerical solution of shock entropy wave interaction CFL = 0.8, = 800 using 

shock switch parameters e = 1 x 10“®,h = 0.8: high resolution of smooth extrema and steep 
gradient region. 


6.4.2 


Sod test tube 


{p,u,p) = 


(1 %/m^, 0 m/s, 100, 000 A^/m^) x<0 

(0.125 kg/m^,0 m/s, 10,000 N/m?) x > 0; 


(76) 


This test problem has no sonic point but the contact and shock are very close which cause a 
smeared approximation to the middle contact discontinuity. In Figure [T^ numerical results are 
given and for different choice of shock switch parameters and compared with FLIC. Results show 
that proposed FLWBW-FLIC crisply captures the smooth rarefaction and contact discontinuity 
and shock more accurately than high order TVD scheme FLIC with Minbee limiter. 


6.4.3 Lax Tube 

{p,u,p) = 


(0.445 kg/m?, 0.698 m/s, 3.528 N/m?) x < 1, 
(0.5 fe(7/m^, 0 m/s, 0.571 A^/m^) x > 1; 


(77) 


Compared to Sod tube the shock in this case is very strong and mostly use to check the robustness 
of any schemes. In Figure 17 numerical results obtained by FLWBW-FLIC are given. It can be 
seen that that the method capture the contact and the rarefaction wave with higher resolution 
shock compared to FLIC for various choices of shock parameters. 
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(a) (b) 

Figure 16: Solution Sod shock tube, N = 100, CFL = 0.5 after 67 time steps at T = 0.01 
second with shock switch parameters e = 1 x 10“®, (a) 6 = 0.4, (b) 5 = 0.9 


6.4.4 Laney Test 


{p,u,p) = 


{1 kg/m^, 0 m/s,100, 000 N/m'^) x<0, 

{0.01 kg/m^, 0 m/s,1,000 N/m‘^) x > 0; 


X G [-10,15]. 


(78) 


In this test the density and pressure state on the right side of initial discontinuity is much 
smaller compared to the left state. Therefore computationally, even small oscillations can lead to 
negative density or pressure which results in to non-physical imaginary speed of sound c = 

This makes it an important test to check the non-oscillatory nature of any numerical scheme. In 
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(a) 


(b) 






X (m) 


(c) 


(d) 






Figure 17: Density and energy plots Lax Shock tube case 2, N = 200, CFL = 0.8 after 187 
time steps at T = 0.32 second with shock switch parameters in row (a) e = 1 x 10“^, (5 = 0.2, 
(b) e = 1 X 10-2, (5 = 0.8, (c) e = 1 x 10“^, 6 = 0.4 and (d) e = 1 x 10-^, 6 = 0.8. 


Figure the density and pressure plots obtained by FLWBW-FLIC are given and compared 
for shock switch parameters e, 5. 

Numerical results for ID shock tube test problems show that the corners of rarefaction and 
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(a) 


(b) 


Figure 18: Solution Laney test, N = 200, C-FL = 0.8 after 111 time steps at T = 0.01 second 
with shock switch parameters e = le — 8 (a) 5 = 0.6 (b) 6 = 0.9 


contact discontinuities are better resolved by FLWBW-FLIC compared to centred flux limiter 
shock capturing TVD scheme FLIC. The shock is also crisply captured compared to to FLIC 
with Minbee limiter. The hybrid scheme presented in section 5.5 is robust and works for different 
choices of shock parameters e, 5. 
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6.5 2D Euler Equation 


Consider the two-dimensional 
system 


Euler equations for compressible gas dynamics defined by the 


Ut + F{U)^ + G{U)y = Q, 


(79) 


where 



(^ \ 


( pu \ 


( pv \ 

u = 

pu 

, F = 

pu^ + p 

, G = 

puv 


pv 


puv 


pv -|- p 


\ e y 


\ u{e+p) y 


\ v{e+p) ) 


(80) 


Here p is the density, u and v are velocity components in x and y direction respectively, p is the 
pressure and e is the energy defined by. 


e = 


p p{v? + v'^) 


7 - 1 




(81) 


The Riemann problem for Euler equation (79) can be defined by considering the constant ini¬ 
tial data in each quadrant of unit square [0,1] x [0,1] with center at (0.5,0.5). More precisely 
consider (79) with initial data 


{p,p,u,v){x,y,t = 0 ) = < 


(82) 


(j)i,pi,ui,vi) if a: > 0.5, y > 0.5 
{P2,P2,U2,V2) if a: < 0.5, y > 0.5 
{P 3 ,P 3 ,U 3 ,V 3 ) if a: < 0.5, y < 0.5 
(p 4 ,/O 4 , ri 4 , U 4 ) if a: > 0.5, y < 0.5 

In the 2D Riemann problems due to complex geometric wave pattern most high resolution 
schemes experience problems in yielding oscillations free crisp resolution to solution profile. Such 
2D Riemann problem are numerically solved by using positive scheme and Riemann solvers free 
central schemes in |21) and [19] respectively. Recently some of the these Riemann problems are 
considered to see the performance of a new finite volume adaptive artificial viscosity method in 
|18) and (with slight changed geometry) a HLL Riemann solver in |43] respectively. In this section 
numerical results are given for twelve configurations. The one dimensional scheme presented 
in section 5.5 is extended to two dimensional Euler equation using the Strang dimension by 
dimension splitting technique. In Figure 19 to Figure 30 the contour plot of density are given 
and compared with FORCE scheme for different test cases. In all the figures contour plot by 
FORCE is given in column (a), and by FLWBW-FORCE in column (b) with shock parameters 
e = 1 X 10“®, 6 = 0.6. It is observed that small oscillations can occur with if a flux limited TVD 
scheme is used in hybrid scheme in section 5.5 (Results are not shown here). Note that it is in 


agreement with the comments in jT^, that a over compressive Minmod type limiter can lead to 
spurious oscillations for the schemes proposed therein. 


Configuration 1 


Pi 

Pi 

Ul 

Vl 


1 

P2 

= 0.4 

P3 

1 

P2 

= 0.5197 

P3 

0 

U2 

= -0.7259 

U3 

0 

V2 

= 0 

V3 


0.0439 p4 = 0.15 
0.1072 p4 = 0.2579 
-0.7259 Ui = 0.0 
-1.4045 va = -1.4045 


(83) 


Configuration 2 


Pi 

= 1 

P2 

= 0.4 

PS = 1 

P4 

= 0.4 

pi 

= 1 

P2 

= 0.5197 

P3 = l 

P4 

= 0.5197 

Ul 

= 0 

U2 

= -0.7259 

U3 = -0.7259 

Ui 

= 0 

Vl 

= 0 

V2 

= 0 

U 3 = -0.7259 

Vi 

= -0.7259 


( 84 ) 
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(a) (b) 

Figure 19: Configuration 1: Density contour plot (30 lines) at T = 0.2, sharp resolution for 
rarefaction and ripple in lower rarefaction are sharply resolved in (b) compared to FORCE (a). 




(a) (b) 

Figure 20: Configuration 2: Density contour plot (30 lines) at T = 0.25 non-oscillatory sharp 
resolution for rarefaction waves and corners in (b) compared to FORCE in (a) 


Configuration 3 


Pi 

= 1.5 

P2 

= 0.3 

P3 

pi 

= 1.5 

p2 

= 0.5323 

P3 

Ul 

= 0 

U2 

= 1.206 

U3 

Vl 

= 0 

V2 

= 0 

V 3 


0.029 p4 = 0.3 
0.138 pA = 0.5323 
1.206 ^4 = 0 
1.206 V4 = 1.206 


Configuration 4 


Pi = 1.1 p2 = 0.35 
= 1.1 p2 = 0.5065 
ui = 0 U2 = 0.8939 

Ul = 0 V2 = 0 


P3 = 1-1 
P3 = 1.1 

Us = 0.8939 
fs = 0.8939 


P4 = 0.35 
P4 = 0.5065 

tt4 = 0 

V4 = 0.8939 


Configuration 5 


pi = l P2 = 1 P3 = 1 P4 = 1 

PI = 1 /02 = 2 P3 = A /54 = 3 

ui = -0.75 U 2 = -0.75 U 3 = 0.75 U 4 = 0.75 


ui = —0.5 V 2 = 0.5 ^3 = 0.5 V 4 = —0.5 


(85) 


( 86 ) 


( 87 ) 
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(a) (b) 


Figure 21: Configuration 3: Density contour plot (32 lines), All four shocks are sharply resolved 
with out induced oscillations in FLWBW-FORCE in (b) compared to FORCE in (a). 



(a) (b) 


Figure 22: Configuration 4: Density contour plot (30 lines) non-oscillatory, FLWBW-FORCE 
yield sharp resolution for shocks in (b) compared to FORCE in (a) 




(a) (b) 

Figure 23: Configuration 5: Density plot (30 Lines) all four contacts are poorly resolved by 
FORCE (a), though FLWBW-FORCE yield sharp resolution to contacts with out oscillations 
in (b). 


Configuration 6 


Pi = 1 P2 = 1 4*3 = 1 4'4 = 1 

Pi = 1 P2 = 2 P3 = 1 P4 = 3 

ui = 0.75 U2 = 0.75 U3 = -0.75 m = -0.75 

ui = —0.5 V2 = 0.5 fs = 0.5 V4 = —0.5 
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(a) (b) 

Figure 24: Configuration 6 : The ripples in NE and SW quadrants are captured with comparable 
resolution with the one in |19) using FLWBW-FORCE (b) though the resolution for contacts is 
little diffusive but much sharper compared to FORCE (a). 


Configuration 7 


Pi = 1 P2 = 0.4 p 3 = 0.4 

Pi = 1 P2 = 0.5197 p 3 = 0.8 

ni = 0.1 U2 = —0.6259 M 3 = 0.1 

Ml = 0.1 V2 = 0.1 M3 = 0.1 


Pi = 0.4 
Pi = 0.5197 

M4 = 0.1 

M 4 = -0.6259 


(89) 



(a) 



(b) 


Figure 25: Configuration 7: The contacts in South and West are crisply resolved by FLWBW- 
FORCE (b). Moreover the rarefaction corners in NE quadrants are significantly sharper than 
FORCE (a). 


Configuration 8 


Pi = 0.4 
Pi = 0.5197 

Ml = 0.1 

Ml = 0.1 


P2 = 1 
P2 = 1 

U2 = -0.6259 

V2 = 0.1 


P3 = 1 

P3 = 0.8 
M3 = 0.1 
M3 = 0.1 


Pi = l 
Pi = l 
M4 = 0.1 

M4 = -0.6259 


(90) 
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(a) (b) 

Figure 26: Configuration 8: The semi-circular wave front in NE is sharply resolved and resolu¬ 
tion for contacts are comparable with the one in |19) . 


Configuration 9 


Pi = 1 P2 = 1 
pl = l p2 = 2 
tti = 0 U2 = 0.0 

vi = 0.3 V2 = —0.3 


P3 = 0.4 
p3 = 1.039 

Us = 0 

ns = -0.8133 


P4 = 0.4 
pA = 0.5197 

U4 = 0.0 

Vi = -0.4259 


(91) 




(a) (b) 

Figure 27: Configuration 9: Again rarefaction and corners are resolved sharper than FORCE 
(a). The vertical contact is crisply captured by FLWBW-FORCE. 


Configuration 10 


Pi = 1 
pi = 1 
ui = 0 

VI = 0.4297 


P2 = 1 
P2 = 0.5 
^2 = 0 
V2 = 0.6076 


Ps = 0.3333 

PS = 0.2281 
Us = 0 

ns = —0.6076 


P4 = 0.3333 
P4 = 0.4562 
n4 = 0 

n4 = -0.4297 


Configuration 11 


Pi = 1 p2 = 0.4 Ps = 0.4 

Pi = 1 P2 = 0.5313 Ps = 0.8 

ni = 0.1 U2 = 0.8276 ns = 0.1 

ni = 0.0 V2 = 0.0 ns = 0.0 


P4 = 0.4 
P4 = 0.5313 
n4 = 0.1 
n4 = 0.7276 


(92) 


(93) 


33 


































































































(a) (b) 

Figure 28: Configuration 10: Resolution of contacts by FLWBW-FORCE is comparable with 
the one in |19| . 




(a) (b) 

Figure 29: Configuration 11: Resolution of contact is SW quadrant, Shocks in SE and NW 
quadrants by FLWBW-FORCE (b) is better than that of the one in |191121] , 


Configuration 12 


Pi = 0.4 
pi = 0.5313 

ui = 0.1 
vi = 0.1 


P2 = l 

p2 = 1.0222 

U2 = —0.6179 

V2 = 0.1 


P3 = 1-0 
p3 = 0.8 

Us = 0.1 
Vs = 0.1 


P4 = 1.0 
p4 = 1.0 
tt4 = 0.1 

Vi = 0.8276 


(94) 


Comments 


7 Conclusion and Future work 

In this work LMP/TVD bounds are obtained for uniformly second order accurate schemes in 
non-conservative form. These bound show that higher than second order TVD accuracy can 
be achieved at extrema and steep gradient region in limiting sense i.e., when r —)> 0“. Based 
on the LMP/TVD bounds hybrid local maximum principle satisfying schemes are constructed 
and applied on various benchmark test problems. Numerical results show improvement in TVD 
approximation of solution region with extreme points, smooth rarefaction as well contact dis¬ 
continuity compared to existing higher order TVD method. For a separate work, the focus is on 
TVD bounds for multi-step methods and efficient use of a shock detector. As, the algorithm |5.4| 
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(a) (b) 

Figure 30: Configuration 12: FLWBW-FORCE recovers the ripples between NE shock and 
contact waves. The resolution for shock and contacts are comparable with the second order 
scheme results in |19) . 


recovered the shock at right location for scalar case, it would be interesting to devise a hybrid 
scheme for system by modifying the wave speed choice in section 5.5, with out a shock detector. 
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